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The actin cytoskeleton of adherent tissue ceUs often condenses into filament bundles con- 
tracted by myosin motors, so-called stress fibers, which play a crucial role in the mechanical 
interaction of cells with their environment. Stress fibers are usually attached to their en- 
vironment at the endpoints, but possibly also along their whole length. We introduce a 
theoretical model for such contractile filament bundles which combines passive viscoelastic- 
ity with active contractility. The model equations are solved analytically for two different 
types of boundary conditions. A free boundary corresponds to stress fiber contraction dy- 
namics after laser surgery and results in good agreement with experimental data. Imposing 
cyclic varying boundary forces allows us to calculate the complex modulus of a single stress 
fiber. 

PACS numbers: 87.10. -he, 87.16.Ln, 87.17.Rt 



INTRODUCTION 



The actin cytoskeleton is a dynamic filament system used by cells to achieve mechanical strength 
and to generate forces. In response to biochemical or mechanical signals, it switches rapidly between 
different morphologies, including isotropic networks and contractile filament bundles. The isotropic 
state of crosslinked passive actin networks has been studied experimentally in great detail, for 
example with microrheology . Similar approaches have been applied to actively contracting 



' 1Ulrich.Schwarz@bioquant.uni-heidelberg.de| 



2 



laser 
cutting 



micromanipulator 




focal 
adhesion 



FIG. 1: (Color online) Contractile filament bundles are very prominent in cell adhesion. Stress fibers typically 
connect two focal adhesions and through their contraction, the cell can probe the mechanical properties of 
the substrate. The viscoelastic properties of stress fibers can be probed by laser cutting or by cyclic loading 
through a micromanipulator. 

actin networks 3, 4] and live cells However, less attention has been paid to the mechanical 

response of the other prominent morphology of the actin cytoskeleton, namely the contractile 
actin bundles, which in mature adhesion appear as so-called stress fibers. During recent years. 



)ut also for 
g. Thus it 



it has become clear that stress fibers play a crucial role not only for cell mechanics, 
the way adherent tissue cells sense the mechanical properties of their environment [7 
is important to understand how passive viscoelasticity and active contractility conspire in stress 
fibers. 

Stress fibers are often mechanically anchored to sites of cell-matrix adhesion, are contracted by 
non-muscle myosin II motors and have a sarcomeric structure similar to muscle 



IC 



Hi, as shown 



schematically in Fig. [TJ However, their detailed molecular structure is much less ordered than in 
muscle. In particular, stress fibers in live cells continuously grow out of the focal adhesions [l^ 
and tend to tear themselves apart under the self-generated stress [isl. Up to now, the mechanical 



response of stress fibers has been measured mainly isolated from cells [14 
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Recently, pulsed 



19l ]. By using the 



lasers have been employed to disrupt single stress fibers in living cells 
intrinsic sarcomeric pattern or an artificial pattern bleached into the fluorescently labeled stress 
fibers, the contraction dynamics of dissected actin stress fibers has been resolved with high spatial 
and temporal resolution along their whole length [1^. These experiments showed that dissected 



stress fibers contract non-uniformly and that the total contraction length saturates for long fibers, 
suggesting that stress fibers in adherent cells are not only attached at their endpoints, but also 
along their whole length. In the same study, cyclic forces have been applied to stress fibers by an 
AFM cantilever, mimicking physiological conditions like in heart, vasculature or gut. Fig. [1] shows 
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schematically how laser cutting and micromanipulation are applied to an adherent cell. 

Early theory work on stress fibers focused on the dynamics of self-assembly leading to a sta- 
ble contractile state 



20, 



2ll ]. Later more detailed mechanical models have been developed and 



parametrized by experimental data 



-27|. Here we investigate a generic continuum model for 
the mechanics of contractile filament bundles and show that it can be solved analytically for the 
boundary conditions corresponding to stress fiber laser nanosurgery and cyclic pulling experiments. 
Our analytical results can be easily used for analyzing experimental data. For relaxation dynamics 
after laser cutting, our model predicts unexpected oscillations. We reevaluate data obtained earlier 
from laser cutting experiments [3] and indeed find evidence for the predicted oscillations. 

This paper is organized as follows. In Sec. |TI]we introduce our continuum model, including the 
central stress fiber equation, Eq. Q. The stress fiber equation is a partial differential equation with 
mixed spatial and temporal derivatives. In order to solve it analytically, in Sec. IIIII we discretize 
this equation in space. This results in a system of ordinary differential equations, which can be 
solved in closed form by an eigenvalue analysis. In Sec. IIVI we take the continuum limit of this 
solution, thus arriving at the general solution of the continuum model. This general solution is 
given in Eq. (j40p . with the corresponding spectrum of retardation times given in Eq. (j38|) . In SeclVl 
and Sec. IVH we specify and discuss the general solution for the boundary conditions appropriate 
for laser cutting and cyclic loading, respectively. In Sec. IVIH we close with a discussion. 

II. MODEL DEFINITION AND SOLUTION 

We model the effectively one-dimensional stress fiber as a viscoelastic material which is subject 
to active myosin contraction forces and which interacts viscoelasticly with its surrounding. In the 
framework of continuum mechanics, the fiber internal viscoelastic stress is given by the viscoelastic 
constitutive equation 



a{t) = f Gint{t-t')e{t')dt' , (1) 

where a = Oxx and e = ^xx = dxU denote the relevant components of the stress and strain 
tensors, respectively. u{x,t) denotes the displacement along the fiber and Gint is the internal stress 
relaxation function. In addition to the viscoelastic stress, the fiber is subject to myosin contractile 
stress am, which we characterize by a linear stress-strain rate relation am = (7^(1 + 5xn/eo)- ^0 
denotes the strain rate of an unloaded fiber and ag is the maximal stress that the molecular motors 
generate under stalling conditions. In addition to the fiber internal stresses, viscoelastic interactions 
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with the surrounding lead to body forces, fext, that act over a characteristic length a along the 
fiber and resist the fiber movement: 

1 /■* 

fext = — Gecct{t-t')u{t')dt' . (2) 

J ~oo 

Here, Gext is the external stress relaxation function. In the following we assume that both internal 
and external stress relaxation functions have the characteristics of a Kelvin- Voigt material: 

G{t) = Ke{t) + 2ri6{t) (3) 

where K and ij are elastic and viscous parameters, respectively, and 6{t) and 5{t) denote the 
Heaviside step and Dirac delta function, respectively. The chosen Kelvin- Voigt model is the sim- 
plest model for a viscoelastic solid that can carry load at constant deformation over a long time. 
Note that Kpp^t represents the elastic foundation of the stress fibers revealed by the laser cutting 



18l |. while r]ext represents dissipative interactions between the moving fiber and the 



experiments 
cytoplasm. 

Our central equation (the stress fiber equation) follows from mechanical equilibrium, 
dx (o" + Cm) + fext = 0, which results in the following partial differential equation: 

dlii + d^u-ru- Ku = . (4) 

This equation has been written in non-dimensional form using the typical length scale a, the time 
scale T = r]int/Kint + crs/{eoKint), the force scale /o = Kint, the non-dimensional ratio of viscosities 
r = arjext/ {Vint + Cs/eo) and the non-dimensional ratio of stiffnesses k = aK^xt/ Kint- Eq. (jlj) 
has been derived before via a different route, namely as the continuum limit of a discrete model 
representing the force balance in each sarcomeric element of a discrete model Q, S^l- However, 
the pure continuum viewpoint taken here seems at least equally valid, because stress fibers are more 
disordered than muscle and because the interactions with the environment represented by -ftText and 
r]ext are expected to be continuous along the stress fiber. In general, the stress fiber equation (jH) 



271]. In this paper, we show that it also 



can be solved numerically with finite element techniques 
can be solved analytically. 

In order to solve the stress fiber equation, we have to impose boundary and initial conditions. 
We impose the boundary conditions that the fiber is firmly attached at its left end at x = 0, and 
is pulled with a certain boundary force fb{t) at its right end at x = L/a = I: 

^x(0, t) = and dxiiih t) + dxu{l, t) + fs = h{t) . (5) 
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The boundary condition for x = I describes the balance of forces at the right end of the fiber where 
f{t) := fs — fb{t) is the difference between the myosin staU force fs = cTs/ fo and the externally 
applied boundary force As initial condition, we simply use ti(x,0) = 0, that is vanishing 

displacement. 

Before we derive the general model solution, we briefly discuss the special case n = V. This 
case can be easily solved and gives first insight into the solution for the displacement field u{x,t). 
With the definition h = u + u, the partial differential equation Eq. @ becomes a homogeneous 
linear ordinary differential equation: 

dlh - Kh = (6) 

with the boundary conditions h{0) = and dxh{l) + f[t) = 0. It can be solved by an exponential 
ansatz and thus leads to a inhomogeneous linear ordinary differential equation u + u = h for u, 
with the initial condition u(x, 0) = 0. The final solution reads 

u{x.t) = - i'^'^lfl, f f{t')e-^'-''Ut' . (7) 



/k cosh(/-^/K) Jq 

Laser cutting experiments correspond to the situation where the externally applied boundary forces 
vanish, that is f{t) = fs = const. Then the integral in Eq. ([7]) is trivial and thus the special case 
K = T leads to a retardation process with a single retardation time r. The largest, always negative 
displacement given by — ^ occurs at x = I, where the fiber was released by the laser cut. The 
magnitude of the displacement decreases exponentially with increasing distance from this point 
and the typical length scale of this decay is given by a/y/K. 

III. SOLUTION OF THE DISCRETIZED MODEL 

In order to find a closed analytical solution for the general stress fiber equation, Eq. we 
discretize our model in space. In order to implement the correct boundary conditions at x = 0, 
it is convenient to symmetrize the system. Thus we consider a doubled model with 2N units 
and 2N + 1 nodes as shown in Fig. [2j Like in the continuum model, internal and external stress 
relaxation are modeled as Kelvin- Voigt-like, that is {hint, ^int) and {k^xt^ lext) are the spring 
stiffness and the viscosity of the internal and external Kelvin- Voigt elements, respectively. Each 
internal Kelvin- Voigt body is also subject to the contractile actomyosin force modeled by a 



linearized force-velocity relationship 



Fm„ = - ^) = + ^{iln - iln-l) ■ (8) 

Vq Vq 
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FIG. 2: Discretizcd model, (a) The filament bmidle is modeled as a linear chain of Kelvin- Voigt bodies, each 
characterized by a spring of stiffness fci„t, a dashpot of viscosity 7i„t and a linear extension a. Actomyosin 
contractility is described by a contractile element with contraction force Fra added to each Kelvin- Voigt 
body in parallel. Viseoelastic interactions between fiber elements and their surrounding are described by an 
additional set of external Kelvin- Voigt bodies with stiffness k^^^t and viscosity 7e2-t. The total fiber length 
is L. Un denotes the displacement of the n-th node, (b) Schematic drawing of the solution for the node 
displacements assuming that both ends are pulled by an external force Fi, > Fg. Since both terminating 
nodes are pulled outward, the solution for the displacements is antisymmetric with respect to the center 
node at ri = 0, which therefore does not move. Thus we obtain the boundary conditions of interest, clamped 
at n = and pulled by _Ff, at n = N. (c) The index n starts counting at the center node. The index j starts 
counting at the node which terminates the fiber at the left. 



Frun is the force exerted by the n-th motor moving with velocity f„. vq is the zero- load or maximum 
motor velocity and is the stall force of the motor. In the final relation we have used that the 
contraction velocity of the n-th motor, Vn, can be related to the rate of elongation of the n-th 
sarcomeric unit as Vn = —{un — Un-i)- 



3C- 



Our model resembles the Kargin-Slonimsky-Rouse (KSR) model for viseoelastic polymers 
33l | , although it is more complicated due to the presence of active stresses and the elastic coupling 
to the environment. The main course of our derivation of the solution for the discrete model 
follows a similar treatment given before for the KSR-model 3J]. The force balance at each node 
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j = 1, . . . , 2N + 1 of the fiber as shown in Fig. [2] reads 

For j = 1: {u2 - ui) - Tiii + {u2 - ui) - kui = -f{t) 

For j = 2, . . . , 2N: [uj+i — 2uj + iij-i) — Tiij + (uj+i — 2uj + uj-i) — kuj = (9) 

For j = 2N + 1: -{U2N+1 - U2n) - ^U2N+1 - iu2N+l - U2n) - KU2N+1 = f{t) . 

We non-dimensionalized time using the time scale r, introduced the non-dimensional parameters 
(k, F) and combine all inhomogeneous boundary terms in the function f{t): 



k. 



■ext 



T 



VOlext 



fit) 



(10) 



Vokint hnt VQ^int + Fs hnt 

It is important to note that Eq. ([9]) is not made non-dimensional in regard to space; this will be 
done later when the continuum limit is performed. 

By taking the difference of subsequent equations in Eq. (j9]) and by introducing the relative 
coordinates yj = Uj^i — Uj, we can write 



M^iscy + M^iasy = fit) 



(11) 



with the 2N x 2N matrix: 







/ 2 + T -1 

-1 2 + T -1 
-1 2 + T 

V ; ; ; 



(12) 



The matrix Meias has the same form as Mmsc, except that k replaces F. In addition we have 
defined the 2A^-dimensional vectors: 

( -fit) \ 



m 



y2 
y2N-i 



and /(t) 







(13) 



We first solve the homogeneous equation. Let A; be an eigenvalue and let vi be the associated 
eigenvector that solves the eigenvalue problem: 



iM^las - \lM^isc)vi = . 



Then the general solution of the homogeneous equation is given by: 

2N 2N 

y(0 = J^QyK*) = J^Q^^ie-^'* 



(14) 



(15) 



=1 



1=1 
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with the eigenvalues and eigenvectors 



K + 4sin2(2(2^) ^ 

A; = K ; and vi 

r + 4sin2(,^^ 



(16) 



It is straight forward to check that Eq. (jl6p is indeed the solution to the eigenvalue problem defined 
by Eq. (jl4p . see the appendix. There we also prove that the 2N eigenvalues are distinct, positive and 



non-zero, and that the eigenvectors are orthogonal and their length is given by vi = \J {2N + l)/2. 
These results validate the form of the homogeneous solution given in Eq. (jlSp . 

In order to determine the solution of the inhomogeneous equation, Eq. (jlip . we use variation 
of the coefficients: 

2N 

y{t) = ^ci{t)vie-^^' . (17) 
1=1 

Inserting this ansatz into the inhomogeneous Eq. (jlip and using the homogeneous solution yields 
2N conditions defining the coefficients ci{t): 



^ci{t)MyiscVie^^'' = fit) 



(18) 



Evaluation of the product M^scVi, rewriting the 2N equations by components, and applying ap- 
propriate addition theorems yields 

27V 



sm( , ) I r + 4sm''(; 



1=1 



'2N + 1' 



'2{2N + 1) 



(19) 



Here the first sinus term is simply the j-th component of the l-th eigenvector. We define a new 
2N X 2N matrix 



2 . , TTlj , 

sml ) 



U 

and a new 2iV-dimensional vector b 

\i 2 ^'^"^ ' ^™ '2(2iV + 1 
With these definitions, Eq. ()19|) can be rewritten as: 



ci{t) r + 4sin2(— — — -) e-^'* 



(20) 



(21) 



Ub{t) = fit) . 



(22) 
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Because U is built up by the normalized and orthogonal eigenvectors, U = I. Moreover it is 
symmetric, thus U = = U^^. Therefore 



b{t) = Uf{t) 



(23) 



The only non-zero components of f{t) are /i = /27V = ~fit)- Therefore the solution for b is given 
by: 

■ ■K2NI \ 

2N + v) 

(24) 

/ 9 / ,.,\ ttI 

-fith 



kit) = -fit) 



2N + 1 



sin( 



2N + 1 



) + sin( 



i+i 



sm 



2N + 1\ ' ' ) ^2N 
We conclude that all even-numbered components of 5 vanish. The coefficients Q(t) are obtained 



by using Eq. (j24|) and integrating Eq. (j2ip : 

( 



Clit) 



fit')e^'''dt' 



2N + lr + 4sm\^^)Jo 
The solution for the relative coordinates follows from Eq. ()17p : 

4 ^ sin(5^)sin(^ ^ 



if I even 



if / odd 



(25) 



yjit) 



2^ + ^ .=i?5,... + 47.nt sin2(^^) Jo 



(26) 



2{2N+1) . 

The actual displacements u^{t) are recovered from the relative coordinates by evaluating the tele- 
scoping sum: 

U2N-\-\-Ux = iu2N+l-U2N) + iu2N - U2N-l) + ••• + (""2 " ^^l) 

V ' 

+ 



y2N 

2N 



y2N-i 



+ 



+ 



yi 



(27) 



Since the solution has to be antisymmetric with respect to the center node at j = -|- 1, compare 
Fig. [21 it must hold true that n27v+i = —ui and more generally n27v+i-fe = ""^i+fci such that for 
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A; = 0, . . . , — 1, the displacements are given by: 

2N--k 



U2N+l-k 2 

j=l+k 

^ 2N-k ^ k 

= 2^ - 2 (28) 

, 2N /2N-k ., k ., 

= - y Q(t)e-^'M y sin(^^)-ysin(^^' 
«=1,3,5,... \ j=i i=i 

In the last step, we used the solution for the relative coordinates given by Eq. (j26p and have 
subsequently reversed the order of summation in both terms. The two sums in parenthesis can be 
further simplified by using the identity 

, , 1 / , , X cosfafn + 1/2)) \ . . 

™0o) = ^ (cot(a/2) ^J^^^ J^ ) . (29) 

Rewriting the result to the index 1 < n < A^, see Fig. [21 we obtain the desired solution of the 
discrete model: 

9 / -,^m-l • 7rn(2m-l) ■ 7r(2m-l) ^ 

u (t) - = V ^"^ ^^+^ ^^+^ / f(t')e~^dt' (30) 

Un{t) - 2/V + l Z^ ■ ^(2m-l) , . ■ 2 7r(2m-l) / ■''^'^''^ '^'^ V'^^'' 

with the retardation times: 

1 r I Icin^ ^(2m-l) 

_ 1 ^ 1 + 4Sm 2(2N+1) , . 

A2m-l,iV K + 4sm^ 2[2N+1) 

Note that Eq. ()30p gives the correct result uq = for the left boundary. For this reason, we can 
extend the range of validity of Eq. (I30p to < n < A^. We also note that the retardation times 
depend on the number of units because the solution describes the movement of a fiber with A'^ 
units which is attached at its left end n = and is pulled at its right end n = N with boundary 
force /fe(t). It is straight forward to confirm the validity of the derived discrete solution, Eq. (I30p 
and Eq. ([3T]) . by inserting it into the discrete model equation, Eq. ([9]). 



IV. CONTINUUM LIMIT OF THE DISCRETIZED MODEL 

The discrete stress fiber model can be transformed to a continuum equation by considering the 
limit A^ — 7- oo while the length L of the fiber is kept constant. In this process, the stress fiber length 
L is subdivided into incremental smaller pieces of length = L/N. Thereby it has to be ensured 
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that the effective viscoelastic properties of the whole fiber are conserved. This is accomphshed 
by re-scahng ah viscoelastic constants in each iteration step with the appropriate scaling factor 



— = ^ according to: 



kN,int = (l>Nhnt and 7Ar,i„t = (pNlint 

(32) 

kN,ext = 4>Nkext and ^N,ext = (p^lext 

To further clarify this procedure consider a single harmonic spring of resting length a and stiffness 
k. This spring is equivalent to two springs of length a/2 and stiffness 2k that are connected in 
series. Here, the scaling factor is 02 = ^ = 2. Thus, the stiffness k^^int in Eq. (j32p represents the 
stiffness of a fiber fragment of length and increases linearly with the number of partitions N , 

ue for the length 



18|. While kN,int 



whereas kint is the reference stiffness of a fiber fragment of length a. A typical val 
scale a would be a = 1 /um, the typical length of sarcomeric units in stress fibers 
increases linearly with the number of partitions, A:7v,ext decreases according to 1/A^. Similarly it 
follows that the viscous parameter ^N,int and jN,ext scale as kN,int and k^^ext-, respectively. The 
non-dimensional parameters and the boundary force scale like: 



-2 



r, KN = (PnI^, fN{t) = (t)Nf{t) . (33) 



i AT = (Pat i , 

We begin the limiting procedure by introducing the continuous spatial variable x = najsf, denoting 
the position of the n-th node within the discrete chain with units. Then Eq. ([9]) yields (also 
compare Fig. [2]) : 

For n = : 
u{0) = 

For n = 1, . . . ,N: 

(34) 

u{x + un) — 2u{x) + tt(x — otv) — T]\fu{x) + u{x + a^v) — 2u{x) + u{x — a^) — k^u = 
For n = N: 

u{L) — u{L — aj\f) + T]yu{L) + u{L) — u{L — a^) + k,nu{L) + /Ar(t) = 

Using the scaling relations for the viscoelastic parameters given in Eq. (I33p and conducting the 
limit iV — 7- oo yields for n = 1, . . . , A^: 

2 1- f u{x + a^) - 2u{x) + u{x - aN)\ r-., . , 

a lim K —Tulx) +... 

N^oo \ a^^ J 

lim ( ^i- + -N)-2u{x)+u{x-a^)\ _ ^^^^^ 

N^oo \ UN ) 



(35) 
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Since oat is a sequence which converges to zero, the hmits define the second derivative of u with 
respect to x. The continuum hmit of the upper equation results in a partial differential equation 
for the displacement u{x,t). The highest order term will contain mixed derivatives in x and t, 
namely, d'^ii. Similarly, the limiting process can be performed for the boundary condition at the 
right end. Note that at this point the spatial variable evaluates to x = Na^ = L: 



a hm r hm ^n(L) + . . . 

a hm ( <L)-u{L-a^)\ ^ ^ a^^^^^ ^ ^^^^ ^ ^ 



(36) 



In each line of the equation, the first limit gives the first derivative of u with respect to x evaluated 
at X = L and the second limit in each line vanishes as a]\f converges to zero. Consequently, in 
the continuum representation, the stresses which originate from shearing the environment cannot 
contribute to the boundary condition. Our continuum model for stress fibers defined by Eq. (jH) 
and Eq. ([5]) is recovered after non-dimensionalizing x, L, u, f using the typical length scale a. 

To obtain a closed solution for the continuous model, we apply the continuum limit to the 
discrete model solution. The limiting procedure is first performed on the retardation times of the 
discrete model given by Eq. ([HT]) : 

Ti^. I ^ -2 7r(2-m-l) p/2 , a at2 ^■„2 7r(2-m-l) 

_ ljV + 4sm 2(2jv+i) _ +4iV sm 2(2jv+i) 
KAT + 4 sm 2(2Ar+l) + 4iV sm 2(2Ar+l) 

Performing the limit N ^ oo yields the retardation times of the continuum model: 

4At/2+7r2(2m-l)2 • ^^^^ 

Since 1 < m < oo, the upper relation defines infinitely many discrete retardation times, non- 
dimensionalized by r. From Eq. (j38p we deduce that the retardation times are bounded by the 
extreme values ri and Tqo according to: 

The first relation holds if k < F, whereas the second holds if k > F. In the special case k = F 
the finite range of possible values collapses to the single retardation time r which leads to a very 
simple form of the analytical solution as shown above with Eq. ([7]). 

The limiting procedure applied in Eq. (I37p can be carried out similarly on the remaining A^- 
dependent terms of the discrete solution given by Eq. (j30p . This yields our central result, i.e. the 
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FIG. 3: First oscillation, (a) If F < k, the displacements of inner fiber segments exhibit damped oscillations 
around their final steady state Uss{x). Here we show a log- log-plot of the time course of the absolute 
difference \u{x, t) — Ugs{x)\ calculated from Eq. (pij) for the position x = 26.4, which is close to the position 
with maximal amplitude and corresponds to the band n = 7 in Fig. [HI The inset gives the amplitudes of 
the first and second oscillation along the fiber, with maxima 236 nm and 9 nm, respectively. The position 
X = 26 A used here is highlighted as dashed line, (b) The difference u{x,t) ~ Uss at the position x — 26.4 is 
shown on a linear scale. Numbering of the extremal values are included for comparison with (a). Parameters 
for (a) and (b) are as in Fig. El (k, /«, F) = (0.028, 0.39, 0) and / = 42.2, a = 1 /im. 



solution for the continuous boundary value problem defined by Eq. Q and Eq. ([5]): 



uix 



(40) 



4r/2 +7r2(2m- 1)2 

Eq. (j40p in combination -with Eq. ()38p is the general solution of our continuum model. We suc- 
cessfuhy checked the validity of our analytical solution by comparision -with a numerical solution 
of Eq. (jH. One big advantage of the analytical solution is that it can be easily used to evaluate 
experimental data. In the following, we will discuss its consequences for the two special cases of 
laser cutting and cyclic loading. 



V. LASER CUTTING 



If a fiber is cut by a train of laser pulses, then there are no external forces acting anymore on 
the free fiber end and fb{t) vanishes. Then Eq. (j40p can be written as: 



(41) 



u{x,t) = ^ Sm{x) (l - e 



m=l 



The solution for the displacement can be understood as a retardation process with infinitely many 
discrete retardation times given by Eq. (|38p and associated, spatially dependent amplitudes 
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Sm{x)- The amplitudes are given by: 



Smix) = 8fsl 




(42) 



The solution for the displacement at x = I, the position of the cut, is particularly simple. At this 
special position, the amplitudes have a linear relation to the corresponding retardation times: 



Since the range of possible retardation times is bounded according to Eq. p9|l . it follows that the 
spectrum at x = I has only negative amplitudes and the resulting solution for the displacement 
at a; = / is always a monotonically decreasing function. However, this is not true for arbitrary x. 
Inspection of Eq. (j42p yields that negative as well as positive amplitudes appear simultaneously. 
Since x = I evaluates the numerator in Eq. (|42p at its maximum, the resulting spectrum constitutes 
a lower bound for the negative amplitudes of the spectra with x ^ I. Similarly, the absolute value 
of Eq. (|43p gives an upper bound for all positive amplitudes. Thus, the retardation spectra with 
X 7^ / oscillate around zero within an envelope for the amplitudes that decays linearly toward zero. 
This can lead to damped oscillations in the displacement of inner fiber bands about their stationary 
value. An representative time course is shown in Fig. [3j The emergence of these oscillations is 
particularly interesting since the stress fiber is modeled in the overdamped limit, that is, inertia 
terms are neglected. We find that these damped oscillations in this inertia-free system occur only 
for T/k < 1, but then for all positions x ^ I. The amplitude of these oscillations reach their 
maximum at distinct positions along the fiber, as shown by the inset to Fig. [3j The location of 
the maxima moves further away from the cut (toward smaller x-values) with increasing order of 
the oscillation. Fig. [3] shows the time course of the displacement aX x = 26.4 which is close to the 
position where the first oscillation reaches its maximum. Using the same parameters as in Fig. [3l 
we show in Fig. [H the time course of the displacement at x = 11.0 where the second oscillation 
reaches its maximum. Since the oscillations are strongly damped, the maximum amplitude of the 
oscillations also decreases with the order. While the maximal amplitude of the first oscillation can 
reach hundreds of nanometers (236 nm at x = 24.9, see inset to Fig. [3]), the maximal amplitude 
of the second oscillation is already much smaller and only of the order of tens of nanometers 
(9nm at X = 11.0, see Fig. U]). Thus, in order to detect the oscillations in experiments, it is 
essential to measure close to where the oscillations reach their respective maximum amplitude. 
We demonstrate this by showing predicted time courses of the difference u(x, t) — Ugs at the two 
positions in Fig. [3] (b) and Fig. [J] (b), respectively. While the first oscillation is most prominent in 



(43) 
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10° I 10^ 



FIG. 4: Second oscillation, (a) Time course of the absolute difference \u{x, t) — Uss \ on a log-log-scale at the 
position X = 11.0, where the amplitude of the second oscillation attains its maximum. The inset again shows 
the maximum amplitude of the first and second oscillation along the fiber, but now the position x = 11.0 
is highlighted as dashed line, (b) The difference u(x,t) — u^s at the position a; = 11.0 is shown on a linear 
scale. Numbering of the extremal values are included for comparison with (a). Parameters used for (a) and 
(b) are the same as in Fig. [3] and Fig.|6l 



Fig. [3] (b), the second oscillations is not detectable at this position. In contrast, the amplitude of 
the first oscillation in Fig. U] (b) is reduced compared to Fig. [3] (b) but the amplitude of the second 
oscillation is much larger and becomes detectable. 

10° 



10-12 
10-2 




-continuum solution 
■ discrete solution 



10- 



10 



10° 



10^ 



FIG. 5: Comparision of t)^Ussix)\ for the continuum model (solid) and for the discrete model (dashed). 
The continuum solution is calculated for a fiber of length / = 42 at position x ~ 26. The discrete solution 
is calculated for a fiber with iV = 42 subunits at node n ~ 26. Both solutions were calculated for the same 
parameters (k, r, /s) = (0.028,0,0.39) as extracted from experimental data. The two solutions agree well 
with each other and both models predict oscillations. 



To show that the oscillations are not an artifact introduced by the continuum limit, we have 
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t[s] 

FIG. 6: (Color online) Experimental results, (a) GFP-actin stress fiber prior to patterning by photo bleach- 
ing, scale bar: 3 fim. (b) Fiber bleached with stripe pattern, (c-e) Stress fiber Is, 30 s and 140 s af- 
ter laser cutting, (f) Time-space kymograph reconstructed from fluorescence intensity profiles along the 
stress fiber. Band positions are extracted by edge-detection (solid lines), scale bars: 30 s and 3 /im. (g) 
Model fit to the displacement data of shown bands (n increasing from bottom to top) with initial positions 
a;n=i,2,5,7 = (42.2, 39.7, 31.8, 26.4)^m yields (k, /,,,r,r) = (0.028, 0.39, 34s, 0.0) with a = LO^m. Note that 
F/k <C 1 and that the experimental data provide evidence for the predicted oscillations, because the curves 
for n=5 and 7 show dips at 8s and 5s after cutting, respectively. 

compared solutions of corresponding continuous and discrete models. Results are shown in Fig. [5] 
where we have used the same parameters as for Fig. [3] and Fig. [H To facilitate comparison 
between continuum and discrete model, solutions are calculated for integer fiber lengths and integer 
positions. We find that continuum and discrete solution agree very well and, most importantly, 
both predict oscillations when T < k. 

Because of the analytical solution Eq. (I4ip . we can easily apply our model to evaluate experi- 
mental data for stress fiber contraction dynamics induced by laser nano-surgery [l^. Briefly, Ptk-2 
cells were tranfected with GFP-actin and a stripe pattern was bleached into their stress fibers. 10 
s later the stress fibers were cut with a laser and their retraction was recorded over several min- 
utes. Kymographs were constructed and for each band, the retraction trace was extracted by edge 
detection. Least-square fitting of the theoretical predictions to four selected bands simultaneously 
was used to estimate the four model parameters (k, /s,r, P). An representative example for the 
outcome of this procedure is shown in Fig. [6] (more examples and the details of our experiments 
are provided in the supplementary material). We find that P = (0.52 it 0.23) • 10~'^ <^ 1 and 
P/k = 0.013 lb 0.021 <C 1 (mean it std, N = 6). This means that in our experiments the second 
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relation of Eq. (j39p applies and that the oscillations demonstrated by Fig. [3] are predicted for this 
experimental system. For the positions corresponding to bands n = 5 and n = 7 our model predicts 
minima at 8 s and 5 s, respectively. Indeed these minima appear as dips in the experimental data 
shown in Fig. [IJg). 



VI. CYCLIC LOADING 



The response of stress fibers to cyclic loading is characterized by the complex modulus which we 
derive from the general solution Eq. (I40p by assuming a cyclic boundary force fb{t) = /s + /oe*^*, 
with a constant offset compensating the stall force of the molecular motors. Evaluation of the 
resulting integral in Eq. (jiO]) yields: 



n(x,t) = 8//oy V4 ^ (e'^'-e~-) . (44) 

Inspection of the time-dependent terms yields that the solution for the displacements approaches 
a harmonic oscillation. The deviations decay exponentially in time, according to e"*/"^™. As a 
consequence, in the limit for large times, the fiber displacements also oscillate with the same 
frequency a; as the force input, but the stationary phase shift between displacements u{x,t) and 
fbit) might vary spatially along the fiber. In the following, we are only interested in the response 
of the fiber as a whole, i.e. we focus on the displacement at x = /. With the above arguments, we 
find in the limit for large times: 

^ J^^4nP + J(2m^J^J^:^ - ' ^^^^ 

The complex modulus, non-dimensionalized by Kint, can be deduced from Eq. (I45p by noting that 
the cyclic force input /qc*'^* and the creep response of the fiber u{l, t) are connected by the inverse 



of the complex modulus 



28|. The expression for the complex modulus can be separated into its 



real and imaginary part, the storage and the loss modulus, respectively: 



V ' ^ V ' (46) 
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FIG. 7: (Color online) Log-log plot of storage modulus (solid) and loss modulus (dashed). Scaling at low 
and high frequencies is shown for both storage and loss modulus. Used parameters are (k, F) = (0.1, 0) and 
I = 30. 



with 



MK/2 + ^2(2m- 1)2 + 1 ' 



(47) 



q[uj) 



81 



UJT„ 



4k/2 + ^2(2m - 1)2 co'^tI + 1 ■ 

An alternative, more concise expression for the complex modulus can be derived by solving 
the Laplace-transformed model equation for the situation of a sudden force application, fbit) = 
fs + fo0{t)j where 0{t) is the unit step function. Solution of this Laplace-transformed boundary 
value problem for u{l, s) = u{l, t)e~^*, with 3 = 7 + iuj, directly yields the Laplace-transformed 
creep compliance, J(s) = u{l,s)/ Jq. It is connected to the complex modulus by: 



G*{(^) = lim 



1 \/l + iuj\/iTuj + K 



y^OsJ{s) tanh ( Z v^^to 



(48) 



Eq. (j46p or Eq. (j48p are equivalent expressions for the complex modulus of the stress fiber model. 
To further study its frequency dependence we use Eq. ()48p . In the special case L/k = 1, it simplifies 
to G*{(^) = {l+iuj)\'^/ tanh(Zy^). The storage modulus becomes a constant, and the loss modulus 
is linearly dependent on the frequency. These are the characteristics of a Kelvin- Voigt body. The 
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more the ratio T / k differs from unity, the larger are the deviations from these simple characteristics. 
To study the general case F/k 7^ 1, consider the limits a; — t- and uj — t- 00. In both limits, the 
stress fiber model again exhibits the characteristics of a Kelvin- Voigt body. The explicit values for 
the limit a; — >• are: 

Q'q = -v/KCOth(/-v/K) 

^ (49) 

Qq{uj) = — csdc? {lyjll) {21k{k — r) + + r) sinh(2/-v/K)) ^■ 

Similarly, in the limit — t- 00, we find: 

g'^= ^ csch2(/^/f ) {21T{T - k) + Vf{T + k) sinh(2/\/f)) 

(50) 

= Vf coth(/Vf) UJ. 

It holds that Q'^o/Gq > 1, with equality for T = k. K similar relation holds for the slope of the loss 
modulus at high and low frequencies. Fig. [7] shows the predicted frequency dependences of Q' and 
Q". 



VII. DISCUSSION 



Here we have presented a complete analytical solution of a generic continuum model for the 
viscoelastic properties of actively contracting filament bundles. Our model contains the most 
important basic features which are known to be involved in the function of stress fibers, namely 
internal viscoelasticity, active contractility by molecular motors, and viscous and elastic coupling 
to the environment. The resulting stress fiber equation, Eq. can be solved with numerical 
methods for partial differential equations. In this paper, we have shown that a general solution 
can be derived by first discretizing the equation in space. In order to implement the correct 
boundary conditions, the system is symmetrized by doubling its size. The resulting system of 
ordinary differential equations leads to an eigenvalue problem which can be solved exactly, leading 
to Eq. ()30p . A continuum limit needs to take care of the appropriate rescaling of the viscoelastic 
parameters and finally leads to the general solution Eq. ()40p for the stress fiber equation. The 
validity of our analytical solution has been successfully checked by comparing it with both the 
discrete and numerical solutions. 

Due to their analytical nature, our results can be easily used to evaluate experimental data. Here 
we have demonstrated this for the case of laser cutting of stress fibers. In an earlier experimental 
study [3], we focused on the movement of the first three bands (n = 1, 2, 3) of the stress fiber. These 



20 



bands are within less than 10 ^im from the fiber tip and thus we did not report the oscillatory feature 
of bands farther away from the cut. After prediction of these oscillations by our analytical results, 
we evaluated the experimental data in this respect and indeed found evidence for their occurance 
(Fig. [6] and supplementary material). This was possible with conventional light microscopy because 
the amplitude of the first oscillation can reach hundreds of nanometers. The amplitude of the second 
oscillation, however, is predicted to be typically on the order of tens of nanometers, which is below 
our resolution limit. In the future, super-resolution microscopy or single particle tracking might 
allow a nanometer-precise validation of our theoretical predictions. 

The extracted parameter values suggest that the frictional coupling between stress fiber and 
cytoplasm, quantified by F, is not relevant in our experiments, and that the retraction dynamics 
is dominated by the elastic foundation quantified by k jisl . However, the elastic coupling to the 
environment might depend on cell type and substrate coating. In fact our findings differ from the 
results of an earlier study, which neglected elastic, but predicted high frictional coupling {25 1. In 
our model, high frictional coupling corresponds to F/k > 1 and thus no oscillations are expected 
in this case. It would be interesting to cut stress fibers in cells grown on micro-patterned surfaces 
that prevent substrate attachment along the fiber. We then would expect not only the transition 
from elastic to viscous coupling, but also the disappearance of the oscillations. 

As a second application of our theoretical results, we suggest to measure the viscoelastic response 



function Gfy) of single stress fibers. This could be done wit^ 



live cells 



18| or on single stress fibers extracted from cells 



1 AFM or similar setups either on 
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16t |. In this case, our model could 



provide a valuable basis for evaluating changes in the viscoelastic properties of stress fibers induced 
by changes in motor regulation, e.g. by calcium concentration or pharmacological compounds. 

In summary, our analytical results of a generic model open up the perspective of quantitatively 
evaluating the physical properties for any kind of contractile filament bundle. In order to apply this 
approach to more complicated cellular or biomimetic systems, it would be interesting to go beyond 
the one-dimensional geometry of bundles and to also consider higher dimensional arrangements of 
contractile elements [l-^, which could be modeled for example by appropriately modified two- and 
three-dimensional networks 35l437l|. 
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IX. APPENDIX: PROOFS FOR EIGENVALUES AND EIGENVECTORS 



In the main text we have used the eigenvalues and eigenvectors given by Eq. (I16p without 
proving that this system indeed solves the eigenvalue problem defined by Eq. (jl4p . Here, we verify 
the solution to the eigenvalue problem and prove the following properties of the eigenvalues and 
eigenvectors: 

(1) The eigenvalues are distinct, positive and non-zero. 



(2) The eigenvectors are orthogonal and their length is given by vi = \J {2N + l)/2 

In order to verify the given eigenvalues and eigenvectors, we first rewrite the matrix -Mg^as ~ 
as: 

( 2B- A -B --A 
-B 2B-A -B ■■■ 
-B 2B-A ■■■ 

V : •••/ 

Where B = 1 — Xi and A = —k + XiT. By using Eq. ()16p one can express A in terms of B: 



M, 



elas 



(51) 



A = 'iB sin^ 



2(2iV + 1) 

Substitution of this relation into Eq. (j51|) and the application of the addition theorem cos 2a 
1 — 2 sin^ a yields: 



(52) 



/ 2 cos 



B 



— 1 2 cos 




-1 

2N+1 
-1 







-1 

2 COS (jNTT 



\ 



\ 



BM, 



(53) 



To prove Eq. (jl4p it has to be shown that the product Mivi vanishes for all / = 1, . . . , 2N. The 
m-th component of the vector which results from this product is given below. It simplifies to zero 
after application of the addition theorem 2 cos a sin /3 = sin(a + /?) + sin(/3 — a): 

'7r(m-l)A { \ ■ ( '^'m-l \ . /vr(m + l)/\ ^ 

- sm I -^^T — — I + 2 cos I I sm ( | - sm ( / | = (54) 



2iV + 1 



2iV + 1 



2A^+ 1 



2iV + l 
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Thus we have shown that the system of eigenvalues and eigenvector Eq. (jl6p indeed solves the 
eigenvalue problem Eq. ()14p . 

Next we show that the eigenvalues are distinct, positive and non-zero. The fact that the 
eigenvalues are positive and non-zero follows directly by inspection of Eq. (jl6p and by noting that 
all viscoelastic constants are positive. It remains to be shown that there are no multiple eigenvalues. 
This can be seen after reformulating the expression for the eigenvalues as: 

A, = 1 + , (55) 

' r 4 sin^ — 

Since 1 < ^ < 2N, it holds for the argument of the sin- function that < 2{2N+i) ^ §■ 
this interval, the sin-function increases monotonically and is single- valued. For this reason, the 
eigenvalues. A;, are also single- valued. The eigenvalues increase monotonically with / if k < T and 
decrease monotonically for increasing / if the opposite inequality holds. 

Next we show that the eigenvectors are orthogonal and their length is given by vi = 
■sj {2N + 1) /2. Consider the matrix of normalized eigenvectors, U, defined in the main text 
Eq. (|20p . By means of this matrix, the statement to be shown can be recapitulated as if^ U = I 

jn — 5k^m- The second relation follows since U is obviously symmetric. In the following we 
will evaluate the square of the matrix U by components: 

2N 

( UU)k,m. = ^ UkJ Uj^rn 

2 X ^ vrfcj TTjm 



— E 



^ TTjik — m) T^jik + m] 

cos — cos 



2N + 2N + 1 2N + 1 

j=i 

The finite sums over the cos-functions can be evaluated by expressing it in terms of exponential 
functions. The used identity is: 

^ 1 , sin(a(n + l/2)) ^^ 

g ^°^(^") = sin(a/2) J (57) 

Application of Eq. (j57p in order to simplify Eq. (j56p finally yields: 

/TTTTN 1 ( ■ //I \ \ ik — m)TT 

[ UU)k,m = , sm {{k - mjvr) cot ■ 



2{2N + 1) V 2(2iV-M) 
-sin((fc + ,n)vr)coti|^ (58) 

-|- cos((A; -1- m)7r) — cos((/c — m)-K)) 
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There are two cases, namely k = m and k ^ m, that have to be considered. 

First assume that k = m. In this case, the last two terms in Eq. (I58p just cancel out each other. 
The sin-function in the second term evaluates to zero while the cotangent gives a finite value: since 
< 2(2N+i) ^ ^' singularities are just spared. Thus, also this term vanishes. It is only the 
first term that gives a contribution. Using rHopital's rule, it evaluates to: 

frTTA 1 r sin((A;-m)7r) 



2(2iV + i)_. 

cos ((k — mlvr) 
= lim -T ^ = 1 

^us 2{2Ar+l) 

This result ensures that all diagonal components of are unity. 

Next assume that k ^ m. In this case the first as well as the second term in Eq. ()58p vanish 
since the sin-function evaluates to zero while the co-tangent yields finite values. The last two terms 
further simplify to: 



1 , -l)fc-"i((_l)2m_^^^ 



(60) 



2(2iV + 1)^ ' 

This result ensures that all off-diagonal components of \fi vanish. The combination of Eq. (|59p 
and Eq. (|60p yields ( UU)k^m = ^k,m which was to be demonstrated. Thus we have shown that all 



eigenvectors are of length vi = \l (2N + l)/2 and form a complete orthogonal basis. 
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